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We show that the changes in the electronic density of states (DOS) in graphene induced by 
localized impurities (single, double and multiple) are significantly different from those caused by 
the long-range Coulomb potential. We focus on gapped graphene; a bound state is present within 
the gap, with a certain amount of spectral weight. As the coupling to the impurity increases, the 
state lowers in energy and approaches the lower continuum valence band. The spectral weight of 
this state does not transform into a resonance state in the valence band, so no unusual screening 
effects related to a redistribution of DOS in the continuum is observed. In terms of the continuous 
Dirac limit, this phenomenon is a consequence of the absence of the "potential bump" at infinity, 
which is present in the potential of the effective Schrodinger equation for graphene with long-range 
Coulomb impurity potential. The states induced by short-range impurity scattering in graphene, 
therefore, have distinctly different properties compared with the long-range potential case. These 
properties closely resemble the case of a short-range single impurity in other bipartite lattices, such 
as square, body-centered cubic, and simple cubic lattices. For these bipartite lattices, there is always 
a localized bound state with energy in the band gap for the entire range of on-site coupling strengths. 
In all cases the energy of these states asymptotically approaches the edge of the valence band as the 
magnitude of the coupling strength increases, but never crosses it. 

PACS numbers: 81.05.Uw, 71.55.-i, 71.23.-k 
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I. INTRODUCTION 

Graphene has been used as an analog case to study 
rclativistic phenomena in heavy atoms because the crit- 
ical charge for the Dirac fermions in graphene is of the 
order of unitji Upon introduction of a single charge 
impurity into an otherwise perfect lattice, a bound state 
appears within the mass gap that separates the original 
Dirac cones 2 ^. As one tunes the coupling strength of 
this charge impurity to larger values, eventually a crit- 
ical charge condition is achieved whereby the energy of 
the bound states passes into the continuum. In gapped 
graphene, this crossover has implications for the screen- 
ing properties of the electron gas»i The wave function 
of such a bound state decays exponentially with the dis- 
tance from the charge. As the coupling strength becomes 
supercritical, the screening effect is sufficiently strong 
that the observed effective charge is reduced by almost 
4 elementary units compared with the unscreened casei 
This phenomenon is similar to what happens in atomic 
physics when the elementary charge is around 170 A The 
shape of the cloud of screening charge closely follows the 
shape of the so called "critical state" just before it merges 
with the continuum of the valence band. A couple of 
questions arise: i) Is there a change in the electronic 
structure that depends on the range of potential? ii) 
What are the properties that define whether the impu- 
rity potential strength is critical or not? 

The paper is organized as follows. In the next section 
we obtain the predictions of the long-wavelength Dirac 



formalism in the case of one or few short-range impu- 
rities. In the third section, after we briefly review the 
Green functions on the lattice pertinent to the honey- 
comb lattice, we obtain analytical results for one and two 
impurities on the graphene lattice. In the fourth section, 
we generalize our results for multiple short-ranged impu- 
rities, and then we close with the conclusion. 



II. DIRAC EQUATION WITH SPHERICAL 
WELL 



The motion of an electron with a fixed full (pseudo-spin 
plus orbital) angular momentum in a circularly symmet- 
ric potential U(r) is: 

(E - U(r) - m)A(r) - (d r + -)B(r) = 0, 

r 

(d r --)A{r) + (E-U{r)+m)B(r) = 0, (1) 
r 

where, for definiteness, U(r) = V9(a — r), where 9{x) is a 
Heaviside step function, a is the radius of the well, and V 
is negative (positive) for a well (barrier). The full wave 
function isi: 



*c(r,< 



J_ f e-tV-W+Aif) 



(2) 



Here j is an eigenvalue of angular momentum J z = L z + 
\<r z , and C can refer to either the A or B sublattice. To 



2 



solve Eq. ([T]), one can express B(r) in terms of A(r): 

, J -A(r)-A'(r) 
is — V + m 

In what follows, we will discuss the case when the en- 
ergy level is at either gap edge; for negative (positive) V 
this is the lower (upper) edge. The solutions inside the 
well are: 



A(r) = dVr 



E + m- V 
E-m-V 



J^ 1/2 (^/(E-V) 2 -m 2 r), 



B(r) = C iy /¥j j+1/2 {y/{E-V) 2 -m 2 r). 
The solutions outside the well arc: 



(4) 



A(r) = Chy/iJ E±H1 K {y / m 2 _ E 2 r) 
V E — m ' 



B(r) 



iC 2 ^/rK :j+1 / 2 {\fr, 



Eh 



(5) 



where J a (x) is a Bessel function of the first kind, K a (x) 
is a modified Bessel function of the second kind chosen to 
satisfy the boundary condition at infinity. The equation 
for the energy levels is: 



^Vj^^^E-VY-m^a) 



J j+1/2 (y/(E-V) 2 -m*a) 
K J+1/2 {Vm 2 - E 2 a) 



(6) 



For a solution to exist, the terms on the left-hand side 
and on the right-hand side must be either pure imaginary 
or pure real numbers. This means, that for — m < E < 0, 
and with V < 0, then we require that (E — V) 2 > m 2 
for a solution. We are primarily interested in negative 
values of V. 

We are interested also in what sequence the levels 
will merge into the continuum, depending on their an- 
gular momentum. One can numerically analyze Eq. ([6]) 
but for the sake of clarity we will analyze the effective 
Schrodingcr equation for A(r), obtained from the system 



-A"(r) + 



,3-0 



-(E- V) 2 + m 2 )A{r) = (7) 



This is a wave equation for a particle in a potential with 

■2 

functional form S-^l — (E — V) 2 + m 2 at zero energy. 
Obviously, for states with a higher value of j 2 —j, the po- 
tential curve is higher, particularly near the origin; states 
with a higher value of j will merge into the continuum for 
larger values of \V\. Another way of looking at the prop- 
erties of the wave function as these solutions merge into 
the continuum is to find asymptotic solutions outside the 
well when E — > — m + 0: 



\A(r)\ 2 /r>. 
\B(r)\ 2 /r: 



C 2 2 E- 



1 



-2\Ai 



i- E - m y^rr? 
C 2 1 



E 2 



-2\Aj 



-E 2 r 



r Vm 2 - E 2 



-E 2 i- 



(8) 
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FIG. 1: Asymptotic behavior of wave function close to the 
band edge for quantum well potential. 



Note that the wave function disappears on A-sites and 
^ 2 B becomes non-normalizable when E = — m. 

As we will see from the following, the Dirac approxima- 
tion for graphene gives correct predictions for electronic 
properties at the critical coupling strength, because the 
energies involved are right at the boundary of the valence 
band, where the linear dispersion is most accurate. 



III. ANALYTICAL RESULTS FOR ONE AND 
TWO IMPURITIES ON THE LATTICE 

A. Lattice Green function 

The Hamiltonian of a free electron in the two- 
dimensional graphene lattice, using the tight binding 
model is, 

H Q = -tJ2(4 b i+s + b l+8 a i) +mJ2(aUi - b\h), (9) 

j,<5 i 

where at is the creation operator of an electron on the 
A-atom site labeled j in the honeycomb lattice, and bj+s 
represents the annihilation of an electron on the neigh- 
boring B-atom site labeled j + S. Here 6 denotes the three 
vectors that connect an A-atom site to its three nearest 
neighboring B-atom sites. The parameters t and m rep- 
resent the nearest neighbor hopping probability and the 
mass differentiating the A and B sublattices, respectively. 
The Hamiltonian in k-space can be written as: 



Hn = 



Pk 

- m 



(10) 
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where we have adopted the standard spinor notation for 
the A and B sublattice components of the wave func- 
tion. Here, (f) k = -t e - lk * a {l + 2 cos(k y V3a/2)e lk * a3 / 2 ), 
where a is the distance between neighboring atoms. The 
eigenvalues are 



e k ,± = ±JP{1+Acl+Ac x c y ) 



(11) 



where c x = cos?>k x a/2 and c y = cos^/3k y a/2. The Green 
functions in k-space can be obtained straightforwardly as 



G (k, iujn) 



G° AA (k,iuj n ) G° AB (k,iuj n ) 
G BA (k,iw n ) G° BB (k,iu n ) 



(12) 



where iu) n = iirT(2n — 1) (with n an integer) are the 
Fermion Matsubara frequencies (T is the temperature), 
and each component is given by 

G° AA (k,icj n ) 
G° BB {k,iw n ) 
G AB {k,iw n ) 

G BA {k,iw n ) = -. (13) 
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c k 
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(ioj n H 




- e 2 

c k 




4>% 




(ioj n H 




- f 2 




4>k 






-m) 2 


c k 



We have added /i, the chemical potential, for complete- 
ness, and the superscript '0' serves to remind us that 
these Green functions are applicable to the clean lattice, 
i.e. without impurity scattering. The lattice Green func- 
tions in real space can be obtained by Fourier transform 
from the above Green's functions: 



where 



G AA (l,j,iuJn) G° AB {l,j,iu n ) 
G BA (l,j,iu n ) G° BB (l,j,iu} n ) 



du 



Air 2 

e iti(i x - j x ) e -it>(i x - j x ) e i2v(l y ~jy ) 



dv (14) 



(uj + /i) — m 2 — t 2 (l + 4 cos 2 v + 4 cos u cost;) 

G°i B ( l J> iu n) = _ J~2 J du J dv ( i5 ) 

e iu(l x - 3m ) e -iv(l x -j x ) el 2v(l y -j y )/ 1 + 2e -i« CQSV \ 
X j o • 

(ui + (i) — m 2 — t 2 (l + 4 cos 2 v + 4 cos m cost;) 

The remaining components are readily obtained through 
the relations G BB (l, j,ico n ) — G AA (l,j,iuimrn — > — m) 
and, for the off-diagonal components, G° BA (l, j, iuj n ) = 
[G AB (l,j,iLu n )}*. Through the analytic continuation, 
iio n — > lo + i(5, we obtain the Green functions slightly 
above the real axis, corresponding to the retarded Green 
function. For the particular case of I = j (an on-site 



Green function) , we can obtain the diagonal components 
analytically. From now on we set t=l, which means that 
all energies are measured in units of the hopping energy. 
The result for G AA is (we set fj, = for simplicity, and 
use the definition E 2 = \lo 2 — m 2 |): 
(i) For < E < 1, lo 2 - m 2 < 0, 



Re[G° AA (l,l,co)} 



(lo + m) 



xF( 



7T 1 



-(E 4 + 12E 2 -6) 



2 2 V (VE 2 + l) VE 2 + 9 
lm[G° AA (l,l,u;)] = 0. 

(ii) For < E < 1, lo 2 - m 2 > 
Re[G AA (l M ] = -t±^l 4 



lm[G° AA (l,l,co)] 



V^E vCE+T) 

[ 2^[3-E] (E+lf h 



2(oo + m) 



V3~~E 



16E 



2 y [3-E] {E + iy 



(iii) for 1< E < 3 
Re[G AA (l,l,u,)} = 

lm[G° AA (l,l,Lo)} = 

(iv) for E > 3 



iu+rn) 1 x F{1 (E + SHE-rf 
2u V£ 2 V 16E 



(co + m) 1_ f( _tt J[3-E](E + lY 



2ir y/E 2' 



16E 



n ,J(E + 3)(E-lf 

{e + 3)(e-i) 3) ' 

lm[G° AA (l,l,uj)] = 0. 

In these expressions we have used F(^,x) = f^ 2 (l — 
xsin 2 6*) _1//2 d9, which is a complete elliptic integral of 
the first kind (also denoted as K [x))& 



B. Analytical formalism 



We start with the equation 

G = (I - G°vy 1 G°, 



(16) 
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where G denotes a matrix where different rows (and 
columns) correspond to different lattice sites (an explicit 
example below will make this clearer, see also Ref. M). 

As an example, we consider the specific case with two 
on-site impurities located at the sites labeled and 1 
(without loss of generality, we number the first atom on 
the A sublattice as 0, and we number the first atom on 
the B sublattice as 1; since the lattice is bi-partite, A- 
atoms are denoted by even numbers., and B-atoms are 
denoted by odd numbers). The I — GqV matrix is then 
written explicitly as 



I - G V 



1 - VG° Q0 
VG% 
VG° 20 
VG% 



-VG° 01 
1 - VG° n 
-VG° 21 
-VG° 31 



(17) 



where now the subscripts refer to the two site indices 
(previously written as arguments in, say, Eq. (fl4"|)). and 
V is the strength of the impurity potential at both sites. 
Then 



G , 



E 

k 



{I-G°V)^G%=C k0 Gl/\ 



(18) 



where C is a cofactor in the matrix (fl7j) . and A is the 
determinant of I— G°V. The factor C k j is (— 1) 4+J times 
the determinant of the original matrix excluding the fc-th 
row and j'-th column. Eq. (fT8|) can be expressed as: 



Gjj — CkjG°j/A 



E 

k<l,kjtj 



^kj^kj 



E 

k>l,k^j 



k]^kj 



(19) 



/A. 



Here the number I is given by the number of sites occu- 
pied by an impurity. 

For j > I (away from the impurities), 

T,k>l,kjtj C ko G % = in Ec l- ©• In this case 
Eq. (|T5|) becomes: 



/A. 
(20) 



G jj — E Gk i G l ? - /A — CjjG° J: j + ^ CkjGlj 

k k<l,kjij 

The cofactor Cjj is equal to A, and therefore 



1 'jj — ' '.;.< 



E Gk 3 G kj 



/A. (21) 



Note that when u = —to, Gq (w + iS) ~ 



C. Single impurity scattering 

When there is only one impurity at any A-atom site, 
the Hamiltonian is given by H = Hq+Hi, where H\ = V. 



The corresponding Green functions corresponding to the 
two Hamiltonians are Go(z) = [z — Ho)^ 1 and G(z) — 
(z — H) . By using the T- matrix expansion, we obtain 
the Green function in the presence of a single impurity: 



r _ r , G i0 VG 0j 

ij " ij + 1 - VG° nn 



(22) 



The local density of states (LDOS) at any position on 
the graphene lattice is defined by the imaginary part of 
Green function: 



P{j, 3, w) 



1 



ImGjj (to + i5), 



(23) 



and we have restored the explicit frequency dependence 
for clarity. 

The local density of states at the impurity site is 



p(0,0,w) 



ilm 



f gi 

\l-V( 



VG° ao {uj+rS) 



. The position of the 



bound state in the gap is determined by the solution of 
the equation 1 — V^Gq (u; + iS) = 0. At the lower band 
edge where oj —m, Gq (^ + iS) ^ lo + m = 0. There- 
fore, inspection of the above equation suggests that no 
solution exists, unless V — > — oo. This observation im- 
plies that for any V the bound state will not merge into 
the lower continuum, i.e. no bound state energy crosses 
the edge at uj = —m. For a single impurity on a S-atom 
site, the same remarks apply for a positive impurity po- 
tential, and the upper band edge plays the role previously 
played by the lower band edge. 

To understand how the bound state approaches the 
continuum band edge, we use the asymptotic expansion 
of the complete elliptic integral of the first kinc-i to get 



G 00 (w + id) 



K(uj + to) In |w + m| 
- VK(u + to) In \uj + m\ 



(24) 



where K = -^=- By expanding the Green function near 
u) = — to wc obtain a pole with spectral weight ao = 
—KiJ\ lnwi, where u>\ = u + to is the solution of 



1 - VKuihxux = 0. 



(25) 



It is clear that as u)\ — > a solution will only occur as 
V — > — oo, and the residue corresponding to that solution 
approaches zero. 

At the other extreme, for a very weak (negative) im- 
purity potential, a similar expansion near u> ~ to gives a 
bound state energy asymptotically approaching the up- 
per band edge (let W2 = to — uj): 



-1 



u>2 ~ exp 



2mK\V[ 



(26) 



The spectral weight approaches zero here as well, as 
ao = 2mKuj2 In 2 W2 , which also goes to zero as the upper- 
band edge is approached. 

To summarize the results of this section, we showed 
that, as the (negative) impurity potential decreases from 



5 



zero towards negative infinity, the frequency of the pole 
migrates from +m (upper band edge) to —m (lower band 
edge). As this occurs, the spectral weight first starts from 
zero, grows to some maximum, and then decreases again 
to zero, as the strength of the potential varies from zero 
to negative infinity. 



IV. TWO OR MORE IMPURITY SCATTERING 
A. Exact solution for two impurities 

We now consider the two-impurity case, with one on 
an A-site, (0, 0), and the second on a i?-site i(i x ,i y ). The 
Hamiltonian is 

H = H + V + V 2 = Hi + H 2 , 

where Hi = Hq + V as in the single impurity case. The 
Green functions G°, G 1 and G correspond to H , Hi and 
H, respectively. The T-matrix for this case is 

f = H 2 + H 2 G 1 H 2 + ... 

Therefore, the Green function becomes 

G^VzGl 



Gjk — G] h 



In fact, for the many-impurity case, the T-matrix 
method can be used in a recursive way, 



where 



G n = G r - 



(Tn)i,i 



71-1 



Vn 

l-V n (G n -i) iti 



To compute the local density of states at site (0,0), we 
need Goo(w) (for simplicity we suppress the iS): 



GooM 
= GooM + 



1-V 2 G\ 1 (oj) 

[1 - V 2 G<>i(")} [1 - VG^u)} VV 2 \G^W 
F(w) 



(uj — U!q) + iS ' 



(27) 



where ujq is the energy of the pole, and F(u>) accounts 
for the remaining (non-singular) frequency dependence. 
The actual pole position is the solution of 

[1 - V 2 G°M] [1 - VG&M] = VV 2 \G° i(ujo)\ 2 , 

(28) 

and the spectral weight is given by F(ojq). 
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FIG. 2: (color online) The behavior of LDOS at the valence 
band edge for two impurities. 



by 



Then the Green function in the case V 2 = V is given 



Gqo(w) 



GUl-VG^ + VlG^l 



(i-vakxi-vGM-vtw 

Keeping the leading term for G^ , we find 
G^w) ~ -2mK hxu), 



a 



Cuj x 
2mK' 



(29) 



(30) 



(31) 



where the quantity C = | Cq X | 2 is finite near cu = —m, and 
we get the result near the bottom of the upper band: 

ao = 2mKui 2 In 2 uj 2 , (32) 
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The definitions of u± , u>2 are the same as in the case of sin- 
gle impurity scattering discussed in the previous section. 
The conclusion is the same: as the strength of the im- 
purity interaction increases, the pole moves towards the 
top of the bottom band, but never crosses it. Instead, 
the residue associated with the pole decreases to zero. 
The key difference with the Coulomb case is that these 
are short range impurities, and this leads to qualitatively 
different behaviour. 



The spatial dependance of the second term in Grr is 
determined by the (Gq R ) 2 , where denotes some (ar- 
bitrary chosen) site in the impurity-occupied area. To 
determine the asymptotic behaviour of Gq R , we use the 
method of stationary phase (see problem 5.2 in£), and 
get 



exp(-RV5E) 



(35) 



B. Long-range asymptotes of the Green functions 
for a large number of impurities 

In this subsection we generalize to some extent the 
results obtained in the previous sections. We find the 
long-range asymptotic behavior of the Green function in 
the case of multiple impurities located inside a finite area 
of the graphene sheet. As a particular case this discus- 
sion includes the circular well, discussed in the Dirac ap- 
proximation, at the beginning of the article in Section 
II. Knowing these Green functions' asymptotes we show 
that the spectral weight of the state near the band edge 
(on the verge of entering into the lower continuum) is zero 
and the screening charge is not significantly reshaped by 
this state.— 

We rewrite Eq. (|2"Tj): 



where SE 



G 



RR 



G° 



RR 



E 



rR^rR 



/A. 



(33) 



Here we introduced the following notation: R is the dis- 
tance from the center of the area in which the impurities 
are confined; we will call this area the "potential well"; 
R corresponds to the site index outside the well; a is the 
radius of the well; r and r' are the distances inside the 
circle of radius a, and r is the index of the site inside the 
well. 

The second term in Eq. (|33l) represents the change in- 
duced by the impurity and is responsible for the spectral 
weight of the bound state. Assuming R 3> a, the fol- 
lowing conclusions can be made about the second term. 
A does not depend on R, but G r R depends on R, and 
C r R ~ G° R , as the determinant G r R contains only one 
row with G°. R . 

Considering the sum over r in J2 r <a G r RG° R we see 

that the r dependency comes only from the phase k ■ 
(i^ — in the exponent of the integrands (|14[ [T5|) . At 

the energies near the band edge ui — ¥ — to, the part k ■ ~? 
can be neglected, as the small fc's produce most of the 
integral value and k ■ ~? does not vary significantly near 
the Dirac point. Therefore, 



, G f rG° r ~ (Gq R ) 2 



(34) 



Thus Grr 



cxp(-2RV8E) 



(SEy/ 2 R 

qualitative agreement with the asymptotic behavior pre- 
dicted by the Dirac equation ([8j . As we can see from the 
standard definition of the Green function 7 : 



Grr (lo) 



y, V> ra (R)V> n (R') 



n 



UJ — U) n 



dc 



</> c (R)v>c(R'r 



(36) 

when oj — > u n the r-dependency of Grr(w) coincides 
with i/j n (R,)ip n (R,)* and we can judge if the state that 
is potentially crossing the band edge into the contin- 
uum is normalizable. In our case as SE the sum 
of Grr over R in the plane diverges, and so does the 
sum ^ n (R)ip n (R)* in the infinite lattice for any finite 
normalizing factor. Hence the state merging into the 
continuum can be called non-normalizable or extended 
and as such has zero spectral weight (in the thermody- 
namic limit). This confirms our conclusion that there is 
no such phenomenon like supercritical screening in the 
case of a localized potential (as opposed to a Coulomb 
potential) in graphene. 



V. CONCLUSIONS 

Our results are in agrement with common intuition 
developed in the physics of shallow states in semicon- 
ductors. In particular it is not only a feature peculiar 
to gapped graphene that properties of the states near 
the band edge are strongly dependent on the long range 
"tail" of the impurity potential^ In general, the effective 
mass approach (mostly determined by long-range prop- 
erties of a system) is in good agreement with exact nu- 
merical methods for the energies near the band edge. It 
is exactly near the band edge where the Coulomb tail 
becomes important while it is negligible in computations 
related to deep levels.— As illustrated in Section II, in the 
continuum limit a potential barrier emerges in the effec- 
tive potential at large distances, due to the squaring of 
the Coulomb potential. Because of the two sublattices, a 
similar 'squaring' occurs when the problem is solved on 
a lattice, and lattice Green functions are utilized. 10 
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